Correlating Basal Gene Expression across Chemical Sensitivity Data to Screen for Novel Synergistic Interactors of HDAC Inhibitors in Pancreatic Carcinoma

Pancreatic ductal adenocarcinoma (PDAC) is one of the most aggressive and lethal malignancies. Development of the chemoresistance in the PDAC is one of the key contributors to the poor survival outcomes and the major reason for urgent development of novel pharmacological approaches in a treatment of PDAC. Systematically tailored combination therapy holds the promise for advancing the treatment of PDAC. However, the number of possible combinations of pharmacological agents is too large to be explored experimentally. In respect to the many epigenetic alterations in PDAC, epigenetic drugs including histone deacetylase inhibitors (HDACi) could be seen as the game changers especially in combined therapy settings. In this work, we explored a possibility of using drug-sensitivity data together with the basal gene expression of pancreatic cell lines to predict combinatorial options available for HDACi. Developed bioinformatics screening protocol for predictions of synergistic drug combinations in PDAC identified the sphingolipid signaling pathway with associated downstream effectors as a promising novel targets for future development of multi-target therapeutics or combined therapy with HDACi. Through the experimental validation, we have characterized novel synergism between HDACi and a Rho-associated protein kinase (ROCK) inhibitor RKI-1447, and between HDACi and a sphingosine 1-phosphate (S1P) receptor agonist fingolimod.


Introduction
Pancreatic ductal adenocarcinoma (PDAC) is considered as one of the most aggressive and lethal malignancies. Early onset dedifferentiation and metastasis significantly humper early diagnosis and treatment options which contributes to the devastatingly poor prognosis of the PDAC [1][2][3]. Although the surgical resection followed by adjuvant chemotherapy offers some prospects of long-term survival, over 80% of PDAC patients are diagnosed in terminal stages when distant metastases have already occurred leaving the chemotherapy as the most frequently applied treatment option [4][5][6]. Additionally, the most of the patients who have received surgery suffer from recurrence within a year, which makes the chemotherapy the mainstay of PDAC therapy [6]. Considering the very modest improvement in 5-years survival rate since 1970s from 3% to 11%, there is an urgent need to develop non-surgical therapeutic approaches for effective treatment of pancreatic cancer [7,8].
Gemcitabine, with or without additional chemotherapeutics, is considered as a first line option in chemotherapy of PDAC. Being one of the most chemoresistant cancers, efficacy of currently available chemotherapeutics for PDAC is seriously compromised by a promptly developing chemoresistence [2,6]. Dense stromal environment together with many genetic and epigenetic alterations contribute to the complex mechanisms of the resistance. Inspired by the many epigenetic alterations in the PDAC, preclinical studies investigating the interference of epigenetic therapeutics with epigenetic mechanisms underlying the PDAC resulted in the promising findings. However, increasing amount of evidence suggests that the epigenetic therapeutics offer measurable benefits in PDAC only in combined therapy settings [5,8]. Histone deacetylase inhibitors (HDACi) represent some of the most promising epigenetic therapeutics for the treatment of PDAC [8]. Currently, there are four HDACi approved by US Food and Drug Administration (FDA) while many candidate drugs are in the clinical trials. Despite the many promising preclinical results, the clinical trials investigating the efficacy of HDACi as a single agent in solid tumors gave rather disappointing results [9]. Therefore, the combination of these potent anticancer agents with other chemotherapeutics is proposed to be one of the most promising approaches in the future development of HDACi-based chemotherapy [8,10].
Despite significant breakthroughs in the cancer research during the last decades, high mortality remains a major issue which highlights the urgent need for novel therapeutic approaches to the PDAC. A combination therapy holds the promise for advancing the treatment of PDAC due to the possibility to reduce the development of drug resistance. However, in experimental settings, the number of possible combinations considering the all approved drugs and drug candidates is too large and the success rates of such experiments are low. In an emerging era of personalized medicine, bioinformatics approaches could offer more rational way to select the optimal combination of drugs for the specific sub-population of patients. In this work, we explored possibility of using drug-sensitivity data together with a basal gene expression data on pancreatic cell lines to predict the combinatorial options available for HDACi. Our results identified the sphingolipid signaling pathway with associated downstream effectors as a promising targets for combined therapy with HDACi or future development of multi-target therapeutics. Through the process of experimental validation of the results, novel synergisms between HDACi and a Rho-associated protein kinase (ROCK) inhibitor RKI-1447, and between HDACi and a sphingosine 1-phosphate (S1P) receptor agonist fingolimod have been characterized.

Results and Discussion
Pharmaco-transcriptomics relationships between basal gene expression and chemical sensitivity of pancreatic carcinoma cell lines (CCLs) were investigated using the systematic correlation analysis of large publically available data sets. Correlating the chemical sensitivity data to the basal gene expression was previously shown to be the valuable method for revealing mechanisms of action of small molecules [11], or in defining the therapy response signatures for limited number of HDACi [12]. In this study, the Genomics of Drug Sensitivity in Cancer (GDSC) [13] and the Cancer Therapeutics Response Portal (CTRP) [11] chemical sensitivity data sets for 864 small molecules across up to 38 pancreatic CCL were correlated to the basal gene expression of 19,221 genes obtained from the Cancer Cell Line Encyclopedia (CCLE) [14]. Pearson correlation coefficients have been calculated between area under the curve (AUC) values and basal expression data (expressed as log2RSEM values) for each transcript. Matrix of Pearson correlation coefficients was further normalized using Fisher's z-transformation to adjust for variations in number of tested CCLs for each of analyzed molecules. Generally speaking, transcripts negatively correlated with the chemical sensitivity data could be attributed to the sensitivity of the pancreatic CCLs towards corresponding small molecule, while the positively correlated transcripts could be attributed to the resistance of the CCLs. In alignment to this, transcripts corresponding to the Tukey's outliers of each small molecule were further extracted and stored as a transcriptomics "fingerprints" of drug sensitivity or resistance. Stored transcriptomics "fingerprints" were further used for predictions of possible novel small molecule combinations. Under assumptions that genes from resistance fingerprints are directly involved in a cancer cell machinery for escaping the drug-induced cell death, while genes from sensitivity fingerprints are directly involved in the mechanisms which sensitize cancer cells on therapy, we hypothesized that level of overlap between these fingerprints could be an indicator of synergism between small molecules and could identify novel pharmacological targets for dual targeting in PDAC. In other words, synergistic effects between two pharmacological agents are hypothesized to occur when first pharmacological agent which is able to sensitize the cancer cells utilizing the same, or similar set of genes as the ones involved in developing of resistance on the second pharmacological agent. The general approach of the study is presented on the Figure 1. To further explore this hypothesis and identify novel potential synergistic counterparts to the HDACi, counting of the transcripts' overlaps was performed on the pool of correlation data where only transcripts with significant correlation were considered (see Materials and Methods). Namely, if the transcripts with significant positive correlation of the first small molecule response had overlapped with the transcripts related to negative correlations of the second small molecule response, these associations were counted and stored. Scores for all possible combinations of HDACi and the rest of small molecules were calculated and the probability of each of combination was evaluated as a count of all possible overlaps between significantly correlated transcripts (named "Final_Syn_Score" in the Supplementary Materials, Supplementary Table S1). Final results included 29,233 possible combinations between HDACi and small molecule interactors (see Supplementary  Materials, Supplementary Table S1).
Using the proposed protocol, several known synergistic interactors of HDACi have been identified which added up to the validity of the approach. For example, the synergistic interaction between DNMT1 inhibitor and pan-HDACi was recently described on MIA PaCa2 cells [15]. Sorafenib is one of the kinase inhibitors (BRAF; FLT3; KDR;RAF1) with preclinical data on synergism with HDACi and this combination is currently under clinical investigations for usage in PDAC [8,16]. Another positive examples include combination of HDAC inhibitors with PI3K inhibitor [17], mTOR inhibitor [18], BET (BRD4 in Table 1) inhibitors [19], and combination with gemcitabine [10] (Table 1). However, it should be noted that not all of the inhibitors of abovementioned pharmacological targets were highly scored in our analysis (for full table of predictions see Supplementary Materials, Supplementary Table S1) which could be attributed to the mostly unknown off-target effects of the small molecules. In order to identify the possible signaling pathways which could be targeted to achieve the synergistic effects with HDACi in the PDAC, a pathway enrichment analysis of annotated targets for predicted small molecule counterparts of HDACi was performed. Enrichment analysis revealed sphingolipid signaling pathway (as defined in KEGG database) [20,21] as one of the top-scored pathways involved in the predicted synergies across analyzed datasets of small molecules ( Figure 1 and Supplementary Materials, Supplementary Table S2). Sphingolipid signaling pathway (as defined through the analyzed KEGG and Wiki Pathway databases [22]) encompass several enzymes involved in the sphingolipid metabolism (Sphingosine Kinase (SPHK), Ceramide Kinase (CERK), Ceramide Synthase (CERS), Alkaline Ceramidase (ACER), Acid Ceramidase (ASAH) etc.) as well as sphingosine-1-phosphate membrane receptors (S1PR1-5) and their downstream Rho/ROCK, PI3K/Akt and MAPK pathways. It is interesting to note that current literature data do not recognize any of the elements of sphingolipid signaling (except downstream pathways activated by diverse signals) as potential novel therapeutic targets for achieving the synergism with HDACi.
One of the central pillars of the sphingolipid signaling is so-called sphingolipid rheostat [23]. Sphingolipid rheostat, as a concept within sphingolipid metabolism, could be seen as a dynamical equilibrium between amounts of pro-apoptotic ceramide and mitogenic and anti-apoptotic sphingosine 1-phosphate (S1P). Sphingolipid rheostat, and particularly equilibrium between ceramide and S1P, have been just recently characterized as one of the critical regulators of the pro-and anti-apoptotic signaling in a metabolically dynamic pancreatic cancers [24,25]. Additionally, Speirs et al. identified SPHK1 as a key driver of the conserved S1P: Ceramide imbalance in the pancreatic cancer subcultures and therefore one of the central controlling hubs of the rheostat machinery [24]. The most recent literature data suggested that sphingolipid signaling pathway could be an important novel therapeutic target to suppress the proliferation across pancreatic tumors made up of heterogeneous cell populations, with recognition of potential of SPHK1 inhibitors as an approach to reverse healthy balance of pro-and anti-apoptotic signaling in pancreatic cancers [24,[26][27][28]. In addition to this, our predictions identified synergism between pan-HDAC inhibitor and SPHK1 inhibitor to be highly scored (in first 5% of the ordered list of predictions) (Table 1 and Supplementary Materials, Supplementary Table S1).
To further corroborate our results on potential synergism occurring between HDACi and sphingolipid signaling pathway in pancreatic carcinoma, expression analysis was performed on the data obtained from pancreatic cancer cohort (The Cancer Genome Atlas (TCGA) database-number of samples 179) and healthy control (Genotype-Tissue Expression (GTEx) database-number of sample 171) ( Figure 2). Besides overexpression of class I HDACs in pancreatic carcinoma patients, results indicated significant elevation of some of the key elements of sphingolipid signaling, including the enzymes involved in sphingolipid rheostat. Namely, expression of SPHK1, but not SPHK2, in patient tissue was significantly increased. Additionally, according to the differential expression analysis (Figure 2), expressions of many components of the sphingolipid signaling pathway were significantly perturbed in a pancreatic cancer tissue, with most of them being upregulated. In summary, the expression analysis suggested significant perturbations of the sphingolipid signaling pathway in the pancreatic carcinoma. Considering the overexpression of HDACs in the pancreatic carcinoma and synergy predictions, these findings further supports exploration of dual targeting of HDACs and sphingolipid signaling as novel approach in treatment of PDAC ( Figure 2).
As a final step of our bioinformatics analysis, interaction network of all of the overlapped transcripts was constructed. Under assumption that these overlapped transcripts could be related to the mechanisms of synergy happening on the protein level, network was created utilizing the STRING database which integrates all known and predicted associations between proteins, including physical interactions and functional associations [29].
Aiming to find underlying mechanism of synergy, network was constructed using overlapped transcripts as nodes and STRING interaction annotations as edges. Betweenness Centrality (BC) measurements were performed on created network to identify the most important communication hubs inside constructed network, e.g., the nodes with the largest number of intersections of the shortest communication paths between other nodes. As a node with the highest BC value emerged p53 tumor suppressor, indicating importance of p53 in the networks mediating the synergistic effects between HDACi and modulators of sphingolipid signaling. Through the complex interactions with HDACs including the post-translational regulation of acetylation status, p53 was found to be an important regulator of HDACi-mediated cancer cell death [30]. Additionally, literature data indicates that pro-apoptotic ceramide accumulation is an important downstream mediator of the p53 response further corroborating involvement of p53 in a mechanism of potential synergy between HDAC inhibition and modulation of sphingolipid signaling [31]. Another identified communication hubs with corresponding BC values are presented in the Table 2. Besides the BC, another classical centrality measures such as Degree and Closeness were considered in identification of influential nodes in the constructed network. Moreover, additional centrality measures unambiguously identified p53 as one of the most influential nodes (Supplementary Materials,  Supplementary Table S3). To evaluate predictive power of proposed approach, synergisms between novel HDACi synthesized in our group with modulators of sphingolipid signaling available in One of the pairs listed among 2.3% of sorted solutions was the potential synergy between Rho-associated protein kinase (ROCK) 1/2 inhibitors and pan-HDACi (Table 1 and Supplementary Materials, Supplementary Table S1). Interestingly, ROCK inhibition, besides its effects on invasion and tumor growth, was previously described as valuable "priming" strategy for chemotherapy of PDAC [33,34]. ROCK is also identified as one of the downstream effectors of the sphingolipid S1P signaling [21,35]. In the experimental analysis we have used ROCK1/2 inhibitor RKI-1447 (IC 50 (ROCK1) = 14.5 nM, and IC 50 (ROCK2) = 6.2 nM). Fingolimod-another compound predicted to be synergistic counterpart of pan-HDACi, was available for experimental evaluation in the laboratory at the time of the study. Fingolimod, the FDA-approved agonist and functional antagonist of S1P receptor subtype 1 (S1PR1) [36], was found to be a potent anticancer agent in some PDAC models which could be at least partially attributed to its moderate inhibitory effects of the class I of HDACs [37,38]. Besides S1P receptors, another elements of sphingolipid signaling pathway were recognized as off-targets of fingolimod, including ceramide synthase 2 (CERS2), S1P lyase, and SPHK1. Altogether, fingolimod appears to be an interesting chemical probe to tackle the potential of synergism between HDACi and sphingolipid signaling modulation [36]. However, it is important to note that fingolimod:HDACi pair was scored among first 13% of predictions indicating rather moderate, or less probable, synergy.
To assess the cytotoxic effects of the drug combinations, Chou-Talalay model [39] was used, which requires drugs to be administered at a fixed dose ratio. MIA PaCa-2 and Panc-1 cells were treated with a combination of synthesized HDACi and RKI-1447, or HDACi and fingolimod in a 4 × 4 dose matrix and measured the resultant cell viabilities by the MTT assay (Supplementary Materials, Supplementary Note S1). Isobologram analysis of the drug combination treatment at low concentrations and concentrations up to the IC 50 values revealed synergisms for a combination of compound 6b and RKI-1447 as well as compound 9b and RKI-1447 on MIA PaCa-2 cells (Table 3).   Synergistic effects of the combination treatment with compound 8b and RKI-1447 on cell viability were observed only at the highest tested concentration (Table 3 and Supplementary Materials, Supplementary Note S1). Considering the results on Panc-1 cell line, isobologram analysis of the drug combination treatment revealed slight synergism for the combination of compound 9b and RKI-1447 only at the lowest tested concentrations (Table 3 and Supplementary Materials, Supplementary Note S1). Similarly to the results on synergism between HDACi and RKI-1447, synergism between HDACi and fingolimod was also observed dominantly in 6b and 9b (Table 4 and Supplementary Materials, Supplementary Note S1).   In alignment with the predictions, synergisms between HDACi and fingolimod were moderate and observed mostly for the highest tested concentrations. Although the antagonisms were detected for some of the tested concentration pairs, these antagonisms are of less concern since the most of them were detected for lower fractions of the affected cell viability (third row of Tables 3 and 4

of each pair).
It is important to note that HDAC inhibitors 6b, 8b and 9b were not part of the initial chemical sensitivity dataset used in bioinformatics analysis. Therefore, experimentally detected synergisms for 6b, 8b and 9b indicate that synergisms predicted by bioinformatics analysis are driven mostly by direct HDAC inhibition rather than off-target effects of specific HDAC inhibitors included in the chemical sensitivity dataset.

Data Preparation
Pharmaco-transcriptomics relationships between basal gene expression and chemical sensitivity of pancreatic carcinoma cell lines (CCLs) were investigated using systematic correlation analysis of large publically available data sets: Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/, accessed on 17 June 2022) [13] chemical sensitivity data sets, Cancer Therapeutics Response Portal (CTRP, https://portals. broadinstitute.org/ctrp.v2.1/, accessed on 17 June 2022) [11] chemical sensitivity data set, Cancer Cell Line Encyclopedia (CCLE, https://depmap.org/portal/download/ 22Q2, accessed on 17 June 2022) [14] basal gene expression data set. Pearson correlation coefficients have been calculated between area under the curve (AUC) values and basal expression data (as log2RSEM (RNA-Seq by Expectation Maximization)) for each transcript across all tested cell lines. Matrix of Pearson correlation coefficients was further normalized using Fisher's z-transformation to adjust for variations in number of tested CCLs for each of analyzed molecules. Transcripts corresponding to the Tukey's outlies (1.5 interquartile range) of each small molecule were further extracted and stored as a transcriptomic "fingerprints" of drug resistance or sensitivity.

Synergy Predictions
Each of the transcripts from sensitivity or resistance signatures was further screened on the significance of positive or negative correlation coefficients across all analyzed small molecules. Only significant correlations (p-value < 0.05) were counted and each iteration was recorded as "Count +/−" or "Count −/+". Total score, titled Final_Syn_Score, was calculated by summation across counts. List of annotated targets for all detected synergisms (Final_Syn_Score > 10 and at least one Tukey's outlier transcript shared between small molecules) was further analyzed using the functional enrichment analysis from g:Profiler server (https://biit.cs.ut.ee/gprofiler, accessed on 5 September 2022) [40]. The Genotype-Tissue Expression (GTEx, https://www.gtexportal.org/home/, accessed on 20 September 2022) [41] database and The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga, accessed on 20 September 2022) [42] were used to analyze the gene expression profiles of HDACs and elements of sphingolipid signaling across pancreatic carcinoma patient data. Data was accessed and analyzed through GEPIA web server (http://gepia.cancer-pku. cn/, accessed on 20 September 2022) [43]. Transcripts corresponding to the overlapped signatures identified across sphingolipid signaling pathway were extracted and analyzed through network analysis using the STRING [29] data in Cytoscape 3.9.1 [44]. To identify major contributors to the communication across networks and mechanism of synergy, BC analysis was performed using CytoNCA tool [45]. , and grown as a monolayer in humidified atmosphere of 95% air and 5% CO 2 at 37 • C. in a 5% CO 2 atmosphere at 370C and in humidified incubator.

Quantitative Analysis of Drug Synergy
To determine the synergistic effects of the drug combinations, we performed an MTT viability assay and the combination index method described by Chou and Talalay [39]. Cytotoxic activity of synthesized HDAC inhibitors, ROCK inhibitor RKI-1447 (Selleck Co. (Shanghai, China) and Fingolimod HCl (Sigma-Aldrich, St. Louis, MO, USA) was assessed on MIA PaCa-2 and Panc-1 cells using MTT assay [46]. MIA PaCa-2 (4 × 103 cells/ well) and PANC-1 (5 × 103 cells/well) were treated with synthesized compounds in five different concentrations (100, 50, 25, 12.5, and 6.25 µM), and each concentration is added in five replicates. After 72 h, 20 µL of MTT solution (3-(4, 5-dimethylthiazol-2-yl)-2, 5dyphenyl tetrazolium bromide) (Sigma-Aldrich, St. Louis, MO, USA) was added to each well. Samples were incubated for 4 h, followed by the addition of 100 µL of 10% SDS and incubated at 37 • C. Absorbance at 570 nm was measured the next day. Cell survival (%) was calculated as an absorbance (570 nm) ratio between treated and control cells multiplied by 100. IC50 was defined as the concentration of the agent that inhibited cell survival by 50% compared to the vehicle control.
To determine the synergistic effects of the drug combinations, we performed the combination index method described by Chou and Talalay [39], using the CalcuSyn software (version 2.0 Biosoft, Cambridge, UK). The combination index (CI) was calculated to assess the benefits of combinational treatment in MiaPaCa-2 and PANC-1 cells. A CI equal to 1 indicates an additive effect, CI = 0.3-0.7 indicates synergism, CI = 0.1-0.3 indicates strong synergism, and CI < 0.1 indicates very strong synergism.

Conclusions
Pancreatic ductal adenocarcinoma (PDAC) is considered as one of the most aggressive and lethal malignancies. Due to the high chemoresistance of PDAC as well as difficult early diagnosis which hampers surgical resection, there is an urgent need to develop novel pharmacological approaches for effective treatment of PDAC. In an emerging era of personalized medicine, systematically tailored combination therapy holds the promise for advancing treatment of PDAC, but number of possible combinations considering all approved drugs and drug candidates in clinical trials is too large to be explored experimentally. In respect to the many epigenetic alterations in PDAC, epigenetic drugs including HDAC inhibitors could be seen as game changers. In this work, we explored possibility of using drug-sensitivity data together with basal gene expression data on pancreatic cell lines to predict the combinatorial options available for HDACi and developed bioinformatics screening protocol for predictions of synergy in PDAC. Our results identified sphingolipid signaling pathway with associated downstream effectors as a promising novel target for future development of multi-target therapeutics or combined therapy with HDACi. The experimental validation of the protocol led to the discovery of the novel HDACi-ROCK inhibitor synergism as well as HDACi-fingolimod synergism. Further experimental evaluation of identified synergism between modulation of sphingolipid signaling and HDAC inhibition will be the major direction of our future work. All predictions made through this study are freely available as a part of Supplementary Materials, Supplementary Table S1.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ph16020294/s1, Table S1: The results of synergy predictions; Table S2: The results of pathway enrichment analysis; Table S3: Additional network centrality analysis; Supplementary Note S1: Detailed results of experimental evaluation of synergisms with concentration annotations.

Conflicts of Interest:
The authors declare no conflict of interest.